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Abstract 

The stabihty of periodic arrays of Mahier-Maslowe or Kelvin-Stuart vortices is 
discussed. We derive with the energy- Casimir stabihty method the nonhnear 
stabihty of this solution in the inviscid case as a function of the solution 



> 

m 

' parameters and of the domain size. We exhibit the maximum size of the 

o ■ 

\ domain for which the vortex street is stable. 

in 

\ By adapting a numerical time-stepping code, we calculate the linear sta- 

O ' bility of the Mallier-Maslowe solution in the presence of viscosity and com- 

pensating forcing. Finally, the results are discussed and compared to a recent 
experiment in fluids performed by Tabeling et al. [Europhysics Letters 3, 
459 (1987)]. Electromagnetically driven counter-rotating vortices are unsta- 



I ble above a critical electric current, and give way to co-rotating vortices. The 

importance of the friction at the bottom of the experimental apparatus is also 
discussed. 
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I. INTRODUCTION 



The problem of vortex dynamics is important for the field of chaotic motion and dy- 
namical systems theory, but the discovery of coherent structures in turbulence has fostered 
the hope that the study of vortices will also lead to a better understanding of turbulent 
flows |[T|. The emergence of coherent flow structures is a well-known feature of quasi- 
geostrophic flows [Q], soap films or two-dimensional turbulence and, because of their 
relevance to large-scale geophysical flows, the dynamics of these structures has attracted at- 
tention during the past two decades. Geophysical fluid flows often appear to be dominated 
by a strong but localized vortical structure that lasts for many circulation times even when 
relatively turbulent flows are impinging upon it. Experimental evidence indicates also that 
the planar free shear layer has an organized two-dimensional structure over a wide range of 
Reynolds numbers . When modelling steady-state configurations of geophysical flows as 
solutions of a dynamical system, it is important to analyze their stability in order to see if 
they can describe physically observable situations. Indeed, in the real world, many random 
forces act on the system and the stationary situations we observe must be stable under these 
perturbations. 

The dynamics of coherent structures in two-dimensional geometry has been studied in 
many different experiments using rotating or stratified fluid (see and references therein), 
in a shallow layer of mercury or of electrolyte subjected to a magnetic field. Here let 
us recall the experiment proposed by Tabeling et al [H . A periodic array of counter-rotating 
vortices is driven by electromagnetic forcing. By passing a current through a cell containing 
a solution of sulfuric acid and an array of permanent magnets of alternating polarity at 
the bottom of the cell, the Lorentz force stirs the fluid, producing the vortices. The two- 
dimensionality of the flow is ensured using a shallow fluid layer. The basic results can be 
summarized as follows: 

(i) At low current corresponding to weak forcing and hence low Reynolds number, the 
flow consists of a linear array of counter-rotating vortices. 
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(ii) This state becomes unstable beyond a critical current. The linear array is now 
composed of nonuniform tilted vortices, alternately large and small. 

(iii) A further increase in the current leads to a state with half the number of co-rotating 
vortices as compared to the initial state. 

This experiment has led to a number of studies concerning chaotic regimes fTO]. Let 



us present another point of view. The first question we want to address is the following: 
what is the connection between the patterns of the Navier-Stokes equation and the exact 
solutions of the Euler equation, where solutions of this type are known to exist ? The second 
question is to determine the stability of such coherent structures in the presence of viscosity 
and forcing. 

We have organized the article in the following way. In section |n|, we review some steady- 
state solutions of the two-dimensional inviscid and incompressible fluid motion. We will also 
present the Mallier-Maslowe vortex street that we will study in the remainder of the article. 
Sections [V| and |V] form the heart of the paper. In section [TTTI , we derive analytically 
explicit sufficient conditions for the nonlinear stability estimates in the inviscid case, using 
Casimirs and convexity properties. In section |V|, we discuss the two-dimensional viscous 
flows and present the numerical method used for studying the linear stability. Section [IV (J 



discusses the results. Section p!VD| sets up the correspondence between the results and the 
experiment. Finally, in section ^ the nonlinear evolution of an unstable Mallier-Maslowe 
solution is presented. The results are then discussed in connection with the experiments of 
Tabeling et al 0. 

II. INVISCID FLOW 

For two-dimensional incompressible fluid motion, one obtains from the Navier-Stokes 
equation, by elimination of the pressure, the equation for the stream function ip: 

+ J( V^V^, ^) = Z/VV + Gex. (1) 

at 
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where J{A, B) = A^By — AyB^ is the usual Poisson bracket and v the kinematic viscosity. 
The last term G^xt is due to possible external forcing. 

In the absence of external forcing and viscosity, the vorticity equation reduces to the 
Euler equation: 



dt 

For steady-state flows, this gives 



J(^,V^V^) (2) 



j(^,VV) = . (3) 

Hence it follows that the vorticity uo = —V'^ip is constant along contours of constant stream 
function ip. The study of planar steady-state flows in an ideal incompressible liquid is 
consequently reduced to solving the following equation: 

where F is an arbitrary function. This is a nonlinear elliptic equation for ip and therefore 
admits a continuous multiplicity of solutions associated with the arbitrariness of F{ip). The 
problem of finding steady states of two-dimensional vortices in an inviscid fluid is then 
equivalent to solving the Poisson equation for the electrostatic potential with the charge 
density self-consistently determined. 

The simplest choice for F is a linear function, which already gives many different patterns. 
Indeed Kolmogorov flows, cellular structures with square or hexagonal cells and even quasi- 



crystal patterns are solutions |11] of the Helmholtz equation V ip = Many other 



solutions have been studied in the literature such as the Lamb dipole [p!2| , p!3[ and the non 



symmetric Chaplygin dipolar solutions 



A possible choice for the function F which has been proposed in the literature p!5|-p!7| is 



This uj — ip sinh- relationship is very important because, using a statistical approach, one can 



show that it characterizes the most probable state of a two-dimensional system of ideal 



point vortices. Published data showing the functional dependence of vorticity on stream 
function in long-lived structures, seen in experiments and simulations, seem qualitatively 
consistent with hyperbolic-sine (as in Eq. ^ or exponential profiles (as in Eq. |^ below) 
which follow from entropy maximization |19[ . These studies have been verified by very long- 
time high resolution numerical studies; however, recent work shows some evidence that this 
result depends on the initial conditions [|l3l. Another recent experimental study led to 



the conclusion that the maximum entropy state is unlikely to be reached since the observed 
final states of flow display characteristics that conflict with the statistical theory. 

In the case of an infinite box, the solution to Eq. (M) introduced by Mallier and 



Maslowe 211 which we will discuss in the remainder of the article is: 



, , / cosh — cos X \ ^ , / ecosx\ ,^ 

ipM = In \ = —2 arctanh (61 

y cosh ey + e cos x J \ cosh ey ) 

which describes a stationary pattern in the form of a street of counter-rotating vortices, 
arranged periodically along the x-axis at intervals of tt. A typical solution is shown in 
figure ^(a) for e = 0.3. The parameter e characterizes the vorticity density: when e = ±1, 
we recover the point vortex solution and when e = we have ip = 0. Thus, as e ranges from 
to 1, the flow represented by Eq. (^ ranges from the fluid at rest, to the flow due to a set 
of point vortices on the x-axis. 



A third choice for the function F of Eq. (^) is that of J. T. Stuart |2^, an exponential 
function: 

If we use the change of variable proposed above, this equation is directly related to the 



well-known Liouville equation for a real scalar field ?/'(?/, t), studied by both Liouville pS 



and by Poincare 



The exact nonlinear solution to Eq. (|^): 

'?/'5 = In (cosh y + e cos x) . (8) 
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is called Kelvin- Stuart's cat's eyes p2[ and is illustrated in Fig. |l](b). This solution can 



be derived in an elegant way with the Hirota method, assuming a vortex spacing of 27r. 
This solution is of interest because the solution corresponds qualitatively to the co-rotating 
vortices seen by Tabeling et al. 0. 

Analytical expressions of the co-rotating or counter-rotating vortex streets are especially 
useful for studying the stability of the experimental fluid flows P| described before. It would 
be of interest to find an analytic function connecting the co-rotating and counter-rotating 
solutions. As written by Mallier and Maslowe [^], if such a function exists, it is likely that 



at the best the general equation will be the following: 

It can be shown that, rather than interpolating between Eq. and Eq. (|^), equation 
reduces to the sinh-Poisson equation (|^) for all nonzero values of A. Therefore a solution of 
equation (H) is 

, / cosh — cos , , . 

\ cosh ey + e cos x J 



where e = (l - Ay/-ABj ' and = | arccosh [(5 - A)/2{-AB)2 

Moreover, it is possible to treat the Liouville equation corresponding to the limiting 
case A = 0, as a singular limit of the sinh-Poisson equation by making the substitution 
1 — = X^e'"^^ and ip = u — P in equation (13) , and then taking /? +oo. By carefully 



following what happens in this highly singular limit, Tracy et al p6[ succeeded in exhibiting 
the Liouville solution as a singular limit of the sinh-Poisson solution. However, to our 
knowledge, the function connecting the two solutions has not been found. At this point, let 
us turn our attention to determining the nonlinear stability of the Mallier-Maslowe solutions. 

III. NONLINEAR STABILITY OF THE COUNTER-ROTATING VORTICES 



We are interested in the stability of the Mallier-Maslowe solution (^) in order to explain 
the experimental results presented by Tabeling et al. To establish explicit sufficient 
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stability conditions for all values of e and to study the nonlinear stability of the counter- 
rotating vortices in a domain D of the plane R^, in the former article |2^, we used the total 
energy on this domain 

Hiuj) =11 -|^|^ dx dy 
JJd 2 

= / -ih^ ih .it ds — [[ -ibW'^ih dx dy 
JdD 2 JJd 2 

= - ipio dx dy (11) 
2 JJd 

where we have used the fact that the velocity, and hence vanishes on the boundary 
of D. Since the fluid is inviscid, this quantity is conserved. More generally, one can also 
show [^, that the functionals C^{uj) = J J ^(to) dx dy, called Casimirs, are also conserved 



for any real- valued function $. 

We define a conserved quantity = H + whose functional derivative is: 

DH^{u).5uj = jj^(i){uj) + ^\u))5u dx dy . (12) 

We wish to choose $ so that DHi^{ujM) = 0, where ojm = — V^^/^m and ipM is defined by 
Eq. (D. We obtain 

-<l>"{u;)=^'M = ^=i=— ■ (13) 

y'4u;2 + (1 _ ^2)2 

leading to (1 — e'^) ^ > — $"(c<j) > 0. We will need to bound — $" away from zero. Eq. 
state that 



u! = — - — sinh(4Arcth (yf(a;, y)) with \g(x,y)\ 



e cosx 



cosh ey 



< e (14) 



so that uj is bounded by \lo\ < Umax = sinh(4Arcth e). The calculation 



4cu2 + (1 - £2)2 < (1 _ cosh(4Arcth e) = — ^ (15) 

1 — 



leads to the improved bounds 



' <$"M<T^fe^ • (16) 



1 - £2 - ^ ^ - 1 + + £4 
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However, the bounds in ([T6|) apply only to < cUmax, whereas we will require such 
bounds to hold over the entire real line. We therefore construct a function $ to coincide 
with $ for \uj\ < tUmax and with 



— + a±uo + j3± 



(17) 



for \uj\ > cUjnax- Tlic coustauts a± and (3± are determined by continuity, so that $ is a 
C^-function. 

With these preparations completed, we are ready to define the nonlinear constant of 
motion: 



D 



-5uj (-V^) ^6uj + ^{lom + Sio) - ^{lom) - ^'{ujm).Suj 



dx dy (18) 



and to use it to establish Liapunov stability estimates. Using the bounds (|T6|), we get, 

( 1 + + g4 ) — ^ -H^M + Suj) + ^ujm) + $ {uJm).Suj < ■ 
We introduce /^mm' minimal eigenvalue of the positive operator (— v^), to obtain 

0< JJ^l6u^{-VT'Suj<h-l\\5u;\\h . 



(19) 



(20) 



Combining (p!9[) and (|20|) , we have 
^ l-£2 N 



1 



(21) 



Now consider an initial value of the perturbation Sujq. Since is a conserved quantity, 

1 



L2 



(22) 



This a priori estimate provides suitable norms bounding the growth of disturbances since 
we have finally 

1 



1 - 

1 + 6^2 + £4 min 



\\5uj\\i2 < 



ll'^^olliz 



(23) 
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The solution is nonlinearly stable if the term in brackets is positive. 

Consider for the domain D a rectangular box, with length 2TrN in x and 2£ in the 
minimal eigenvalue of the operator (— V^) is fc^m — + since the eigenfunctions 

vanishing on the boundary are f{x,y) = cos{x/N) sin {-Ky/i). Therefore, we have derived 
a maximum transverse size of the domain D for which the Mallier-Maslowe vortex street is 
nonlinearly stable. The sufficient conditional stability is the following: 



71 1 + + 1 

Figure ^ presents the region of sufficient stability of the counter-rotating vortices in the {£, e) 
plane for = 1. 

IV. VISCOUS AND FORCED FLOW 
A. Introduction 

It would be interesting to extend the previously presented results in the presence of 
viscosity and forcing. Moreover the viscosity imposes a minimum scale (z// sup ci;)^, the 
diffusion length for the eddy turn-over period at the maximum realised vorticity. Thus, a 
small viscosity avoids some difficulties concerning the continuum limit of Euler flow. To 
fully understand the nonlinear evolution, we can follow the time evolution of the system 
from various initial conditions; however, it is useful to obtain the eigenspectrum of the 
steady states, since they are associated with transitions and loss of stability. 

In the viscous case, we have the full two-dimensional Navier-Stokes equation 

^ + J(V2^, ^) = z/VV + (25) 

In what follows, we will choose the external forcing to counterbalance the viscosity: 

Gext = -^^VVm (26) 

for which we shall give a partial justification in section [IV D| . With the choice if ipM 
is a stationary solution to the Euler equation (^), then ipM is also a solution to the the full 
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Navier-Stokes equation (^). Viscosity plays an important role, however, in determining 
the stability of the solution tpM- Our strategy is to linearize the Navier-Stokes equations 
around the steady states and to seek eigenmodes of the linearized equations. The linearized 
equation governing the evolution of a perturbation (p is 

= J{iPm, + J(0, V'V^Af) + //VV • (27) 



B. Numerical procedure 

All of our calculations are performed on a two-dimensional (x, y) plane. In the periodic 
x-direction we use a Fourier representation with A^^^ modes (from 16 to 64); We map y G 
(— cxD, +oo) to (—1, +1) via a tanh mapping with A^j, gridpoints (from 65 to 123). Boundary 
conditions are automatically satisfied in this representation: 0(x + 27r, y) = (p{x, y) and 
d(f)/dy{x,y = ±oo) = 0. 

For stability, the viscous term in Eq. (^) is integrated implicitly by the backward 
Euler scheme. The remaining terms are integrated explicitly. We have 

V'<Pn+l = (/ - uAtV^)-' [V20„ + At (j(V^M, VVn) + J{<Pn, V^^m))] • (28) 

where At is the time step. 

The linear stability of ip^j is determined by the leading eigenvalues (those with greatest 
real part) of the operator on the right hand side of (pTf). The leading eigenvalues of the 
operator of the differential equation (p7|) become the dominant ones (those with largest 
magnitude) of the iterative scheme of Eq. (|28D; fortunately, dominant eigenvalues are those 
most readily calculated by iterative methods. Effectively, exact solution of Eq. (0) would 
require exponentiating the operator on its right hand side, and the numerical method ( pSf ) 
carries out an approximate exponential. 

This can be abbreviated as 

= A ^ VV(i) = e^* V'0(O) = e^"^* V'0(O) = fi" V'0(O) (29) 
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where t = nAt and B = e^^* is approximated by the operator on the right-hand-side of 
Eq. @. 

The block power, or Arnoldi's, method is used in order to find the k leading eigenvalues, 
including complex or multiple eigenvalues, simultaneously, as described by Mamun et Tuck- 
erman [^| and references therein. We first integrate Eq. ( [27| ) for some fairly long period 
of time T in order to purge the vector of the strongly damped eigenmodes which are not 
important for the linear stability study. We then take k additional time steps, creating 
Ml = u{T), . . . ,Uk+i = u{T + kdt). The vectors are orthonormalized, forming a basis for 
what is called the Krylov space. A {k by k) matrix if, which represents the action of B 
on the Krylov space, is generated and diagonalized, yielding eigenvalues and eigenvectors of 
the linear stability problem. The eigenvalues A of A are recovered from those of B (or H) 
by taking their logarithm and dividing by At. 



C. Results of the linear stability analysis 

Let us now study the linear stability of the Mallier-Maslowe vortices i.e., the stability 
of the fiow ipM defined by Eq. (P) with corresponding forcing (p6l). Figure |^ presents the 
real part of the two first eigenvalues obtained with the above method as a function of the 
kinematic viscosity when the vorticity parameter e is 0.3 in Eq. @. 

In the high viscosity regime {v ^ 1), the Jacobian terms can be neglected and equa- 
tion (^) becomes 

^ = ^V'4> (30) 

i.e., a heat equation for the vorticity at z/ ^ +oo. In this limit, the equation is 
independent of the vorticity density parameter e, so the eigenvalues are also independent 
of e. The numerical results confirms that, at sufficiently large viscosity and with the forcing 
chosen according to Eq. (p6D, the Mallier-Maslowe solutions (|^) are stable, since the growth 
rate of perturbations is negative. Around the value of z/ = 0.5, the fiow becomes unstable. 
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The growth rate increases as the viscosity decreases. The numerical method presented in 
section [TV ij| is feasible only for high to moderate viscosities since for low viscosities, stability 
of the explicit part of the numerical scheme requires a very small time step, leading to a 
time-consuming code. However, in the zero viscosity limit, we showed in section |IT1| that the 
instability increases with e. 

In order to more clearly understand the evolution of the most unstable eigenmodes as 
a function of the viscosity, we will study the two particular cases depicted by filled squares 
in Fig. H, one stable {y = 10) and the other unstable (z/ = 0.01). Figure ^(a) depicts the 
least stable eigenvector for e = 0.3 and u = 10. As we see, the eigenvector is reflection- 
symmetric in y and independent of x. This mode resembles a shear layer. The evolution of 
the eigenvalue shown in Fig. ^ attests that the x-independent mode is not the least unstable 
mode for low viscosity; however, it is important for the following discussion to notice that 
this mode is marginally stable for low values of u. Figure ^(b) presents the contour plot 
of the most unstable eigenvector for e = 0.3 and u = 0.01. The eigenvector has the shift- 
and-refiect symmetry 0(x + tt, — y) = y) and has the same periodicity in x as ipM- The 
growth rates of the two different modes cross at around the value u = 0.5, as can be seen in 



With the use of a simplified heuristic model it is possible to understand the two 
modes. Since the flow is mainly present in a confined region, then, for the sake of sim- 
plicity, let us consider the main flow to consist of counter-rotating vortices in a finite 
box: ipQ = A sin{kx) cos{£y) . We approximate the marginally unstable shearing mode by 
01 = Bcos{iy), as suggested by Fig. |(a). With this ansatz, one can then show that the 
Jacobian term in Eq. p7| ) will give 



The interaction of the basic flow ipo with the marginally unstable mode 0i thus gener- 
ates a third term 02 = C cos{kx) sin{2iy), completing the triad. We therefore consider an 
unstable perturbative mode of the form (0i + 02), which will have the same pattern as the 



Fig. |. 




(31) 
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most unstable mode found for v = 0.01 and shown in Fig^(b). 

We can go further and explain the occurrence of the instability. We continue to approx- 
imate the base flow by ipo = Asm{kx) cos{iy), with A flxed, and the perturbation by 

01 + 02 = B{t) cos{iy) + C{t) cos{kx) sm{2iy) . (32) 

We substitute these approximations into Eq. 

aV2(0i + 02) 



dt 

obtaining 



J(V^O, V'(01 + 02)) + J((01 + 02), V>o) + i^V"(01 + 02) , (33) 



- Bf cos{iy) - (A;2 + 4f)C cos{kx) sm{2iy) = — cos(A;x) sm{2iy) 

3 

- -kfA C [cos{£y) + 3 cos(2A;x) cos{iy) + 3 cos{3iy) + cos(2A;x) cos{3iy)] 
+ u[e^ B cos{iy) + (A;2 + Ai^)^ C cos{kx) sin(2£?/)) . (34) 

Projecting onto cos{iy) and cos(A;x) sm{2iy) gives the following Galerkin system for the 
time dependent amplitudes B and C: 

3k 

B= -uf B - C (35a) 
C= -uie + Af) C- ^^f_^^^,^ AB (35b) 

Finally, looking for solutions B = Bq e*** and C = Cq e***, one gets the equation 

+ su{k^ + 5f) + uH^P + Af) - ^^^l^^ = . (36) 

For low values of A (i.e. low value of the intensity of the electric current: see section |1V D| ), 
the initial flow is stable since all solutions of equation (|36|) are negative. In contrast, above 
the threshold value 

= p Vs ' ^^^^ 

one solution of Eq. ( P5| ) is real and positive: we get a stationary bifurcation giving rise 
to an instability of the perturbation (0i + 02) whose pattern coincides with that shown in 
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Fig. 0(b). This simple approach therefore gives a good quahtative understanding of the 
relationship between the marginal and unstable modes, and of the onset of instability. 

Figure ^ shows the evolution of the greatest eigenvalue versus viscosity for three different 
values of £ on a logarithmic scale. We see that the results are independent of e in the 
diffusive regime (high viscosity), as explained in the preceding section. The inset allows us 
to ascertain that the value of A is a linear decreasing function of the viscosity. In the low 
viscosity regime the evolution is qualitatively the same but the curves are distinct. The 
bigger the parameter e, the bigger the leading eigenvalue and, therefore the more unstable 
the flow. One can also verify that the Mallier-Maslowe vortices with e = 0.5 become unstable 
at a critical viscosity which is slightly higher than that corresponding to e = 0.2: the 
critical viscosity is an increasing function of e. 



D. Relation to the experiment 

In the experiment by Tabeling et al. 0], the typical velocity V of the basic regime 
can be found by balancing the forcing with the viscous term, as we have done via our 
assumption ([26|). In dimensional terms, this leads to the relation 

V = ^ (38) 



in which B is the maximum value of the magnetic field, / the intensity of the electric current, 
h the depth of the fluid layer, i the width of the magnet, u the viscosity and p the density 
of the fluid. Because the typical velocity varies linearly with /, it is reasonable to suppose 
that the streamfunction and its derivatives will also increase linearly with I. We thus take 
our streamfunction ip = ^ipM, where ipM is the Malher-Maslowe solution (H) and ^ is a scalar 
which increases with /. 

One can easily check that if ipe is any solution to the stationary solution Euler equa- 
tion (0), then C.'ipe is also a solution. Let us show that the linear stability of all the solutions 
C,ipe is determined by a linear stability analysis of ipe as a function of viscosity. The equation 
governing the evolution of an infinitesimal perturbation to the new inviscid solution is 
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^ = J(e^e, V20) + J(0, + i^VV • (39) 

Dividing by ^, we get 

-^-^ = Ji^Pe, VV) + ^(0, VVe) + ■ (40) 

Thus, studying the hnear stabihty of the solution ^ip^ for viscosity z/ is equivalent to studying 
the linear stability of the solution ipe for viscosity {v/C}) except that the eigenvalue will also 
be modified by ^. 

Using this insight, it is then possible to understand the appearance of the instability. 
Recall from Fig. |^ that the counter-rotating vortex flow ipM is stable for sufficiently high v. 
Thus ^ipM is stable for sufficiently high z//^, i.e., for sufficiently low electric current. In- 
creasing the electric current I in the experiment corresponds to increasing ^ and, therefore, 
to decreasing the "renormalized" viscosity {v/C)- The solution will therefore remain stable 
until the renormalized viscosity reaches the critical viscosity z/^ shown in Fig. ^. Above this 
threshold, the counter-rotating vortices will be unstable and will evolve as presented in the 
next section. The appearance of the instability of the counter-rotating vortices for a high 
enough electric current I is thus explained. 

V. TIME INTEGRATION AND CONCLUSION 

The transitions resulting from the linear instability of the Mallier-Maslowe vortices are 
studied by time-integrating the Navier-Stokes equation (pS]). This numerical experiment is 
constructed to resemble that of Tabeling et al, except that the size of the box is infinite. 
As an initial condition, we add to the Mallier-Maslowe vortex flow ipM a small perturbation 
of the form (l){x,y,t = 0) = exp(— y^) coss, to accelerate the appearance of the possible 



unstable modes. The spatial representation is, as discussed in section [IV B| , with a resolution 
of = 32, Ny = 69 and t/max = 15. The time stepping is carried out according to Eq. (p8|), 
with At = 0.01. 
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If the simulation is carried out with e = 0.3 and z/ = 5, the patterns are stable, confirming 
the linear stability analysis presented in the previous section. When the parameters e and u 
are fixed at 0.3 and 0.01 respectively, the evolution, depicted in Fig. (^, is clearly different. 
We see in Fig. ^(b) that at t = 65 we have a linear array of tilted vortices of positive sign 
and the size of the vortices has doubled as occurred in the experiment. The negative vortices 
have been ejected away from the center of the box. 

If we plot the deviation = {ip — ipM) from the initial condition tpM-, we find in the initial 
stage (at t = 50 for example) the streamfunction presented in figure 0(a), confirming the 
linear stability analysis (see Fig. §(b)). After this linear transient growth, the flow continues 
to evolve until it approaches the pattern depicted in figure |](c). However, the amplitude 
continues to increase with time (see Fig. while preserving the pattern. Figure ^(b) 
depicts the deviation at t = 200 and we note that the pattern around the x-axis resembles 
somewhat the Kelvin-Stuart vortices (see Fig. |l](b)). This explains why the system does not 
reach a final equilibrium state: Since (j) is then itself a solution of the stationary inviscid 
Eq. (|^) J(0, V^0) = 0, the nonlinear saturation effect disappears and continues to grow 
in time. A similar case was found in zero Prandtl number convection, where the linearly 
unstable roll modes are exact nonlinear solutions [PT|. 



One possible reason that our simulation, unlike the experiment of Tabeling et al |^ , does 
not reach a final equilibrium state could be that we considered the flow to be perfectly two- 
dimensional. In the experiment, two- dimensionality is enforced by using a shallow fluid layer: 
a frictional force proportional to the velocity could capture this bottom- friction effect ||32| , |33| . 
The addition of a term —fiV'^ip proportional to the velocity on the right hand side of the 
equation (^5]) will change the evolution of the flow. The eigenvalues in the presence of linear 
friction differ from those of the problem without friction only by a shift (— /i) of the growth 
rate: thus the determination of the dependence of the eigenvalue spectrum on fi does not 
require additional numerical studies and, in addition, the linear friction is always stabilizing. 

We therefore time-integrate this system, including the linear friction term as well as the 
ordinary viscosity. The external force is now chosen as 
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Gext = -z^VVm + /XV'^/'M (41) 

so that, as in the previous section, the viscosity only acts on the perturbation not the basic 
flow. The resulting streamfunction at t = 300 is shown in Fig. ^ for a small value of /i fixed 
at 0.01. Contrary to the case without linear friction (/i = 0, see Fig. |^), the flow attains 
an equilibrium state with co-rotating vortices along the x-axis as demonstrated by the time 
series in Fig. Thus, the linear friction term stabilizes the row of co-rotating vortices as 
was obtained in the experiment. The necessity of this linear friction term in reproducing 
the final state of the experiment could be a reason why the final maximum entropy state is 



not often reached by Marteau and al pO[ in their experiment: in their small 2D lattice of 
electromagnetically forced vortices, the bottom- friction effect should also be important. 

The purpose of this work was to understand and explain the behavior of the instructive 
experiment of Tabeling et al. 0. First, we derived explicitly the nonlinear stability con- 
dition for the counter-rotating vortex solutions in a rectangular box without compensating 
viscosity; the finite-size effects of the box were studied. Then, introducing viscosity and 
compensating forcing, we derived a model for the appearance of the instability when the 
electric current is increased: the renormalized viscosity decreases until it reaches the criti- 
cal viscosity z/^ determined by a numerical linear stability analysis. We found Uc to be an 
increasing function of the parameter e characterizing the vorticity which could explain the 
fact low vorticity- flows (i.e. low e in Eq. (^) are more visible in 2D hydrodynamic flows, 
since they are more stable than high vorticity-flows. Above this threshold, the evolution of 
the system leads to a final equilibrium state similar to that in the experiment if a linear 
term is added to the standard Navier-Stokes equation to reproduce the effects of the friction 
at the bottom of the experimental apparatus. 
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FIGURES 

FIG. 1. Steady Flows. Figure (a) represents the streamlines of the Mallier-Maslowe solution @ 
for e = 0.3. Figure (b) represents the streamlines of the Kelvin-Stuart solution (^) for e = 0.3. 
The dashed curves are negative contour lines and the solid curves are positives ones. 

FIG. 2. The region A defines the domain of sufficient stability of a pair of counter-rotating 
vortices (A^ = 1) in the plane e). i is the transverse size of the box and e characterizes the 
vorticity density. The solid line is defined by equation (^). 

FIG. 3. Eigenvalues. Dependence of the growth rate for the first mode on the viscosity u for the 
counter-rotating vortices when e = 0.3. The squares and the dashed curve correspond to the most 
unstable (or least stable) eigenvalue and the asterisks and the solid curve to the x-independent 
mode. 

FIG. 4. Eigenvectors. The streamlines of the most unstable eigenvector associated with the 
Mallier-Maslowe vortices for e = 0.3 are presented as a surface plot for = 10 in figure (a), 
while figure (b) depicts the contour-plot when v = 0.01. The eigenvalues associated with these 
eigenvectors are A = —0.0115 and A = 0.099, respectively. 

FIG. 5. Dependence of the growth rate on viscosity for three values of e for the counter-rotating 
vortices. The diamonds and the dash-dotted curve correspond to e = 0.2, the squares and the 
dash-triple-dotted curve correspond to e = 0.3 and the triangles and the dashed curve correspond 
to e = 0.5. The asterisks and the solid curve represent the j;-independent mode. A logarithmic 
scale is used for the viscosity; the inset uses a linear scale. 

FIG. 6. Contour plot of the stream function for e = 0.3 and v = 0.01. (a) initial condition, the 
Malher-Maslowe Va/; (b) t = 65; (c) t = 200. 

FIG. 7. Contour plot of the deviation (j) = ip — ipM from the Mallier-Maslowe solution for 
£ = 0.3, u = 0.01 and /i = 0. (a) t = 50; (b) t = 3000. 
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FIG. 8. Evolution of the norm of the deviation ||(^|| = HV' — V'm|| from the MaUier-Maslowe 
solution for e = 0.3, u = 0.01 versus time. The solid curve corresponds to /x = and the dashed 
curve to n = 0.01. 

FIG. 9. Contour plot of the stream function ip ior e = 0.3, v = 0.01 and /x = 0.01 at t = 5000. 
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